Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 1 February 2008 (MN I*TeX style file vl.4) 

CDM-variant cosmological models — I: Simulations and 
preliminary comparisons 



00 



Michael A. K. Gross/'^"*" Rachel S. Somerville/'^'*' Joel R. Primack/'*' 
Jon Holtzman^* and Anatoly Klypin^* 

^Physics Department, University of California, Santa Cruz, CA 95064 USA 

^Earth & Space Data Computing Division, Code 931, NASA/Goddard Space Flight Center, Greenbelt, MD 20771 USA 

^ Racah Institute of Physics, The Hebrew University, Jerusalem 91904, Israel 

^Department of Astronomy, New Mexico State University, Las Cruces, NM 88001 USA 



>. 

cd 



1 February 2008 



(N 



(N 
> 

(N 

(N 

o 



'i 



ABSTRACT 

We present two matched sets of five dissipationless simulations each, including four 
presently favored minimal modifications to the standard cold dark matter (CDM) sce- 
nario. One simulation suite, with a linear box size of 75 h~^ Mpc, is designed for high 
resolution and good statistics on the group/poor cluster scale, and the other, with a 
box size of 300 h~^ Mpc, is designed for good rich cluster statistics. All runs had 57 
million cold particles, and models with massive neutrinos (CHDM-2:/) had an addi- 
tional 113 million hot particles. We consider separately models with massive neutrinos, 
tilt, curvature, and a nonzero cosmological constant (A = iH^VltC} in addition to the 
standard CDM model. We find that the dark matter in each of our tilted fio + ^A = 1 
(TACDM) model with 0.^ = 0.4, our tilted f^o = 1 model (TCDM), and our open 
A = (OCDM) model with Oq — 0.5 has too much small-scale power by a factor of 
^ 2, while CHDM-2z/ and SCDM are acceptable fits. In addition, we take advantage of 
the large dynamic range in detectable halo masses afforded by the combination of the 
two sets of simulations to test the Press-Schechter approximation. We find good fits 
at cluster masses for (5c, g — 1.27-1.35 for a Gaussian filter and (5c, t — 1.57-1.73 for a 
tophat filter. But, when we adjust 5c to obtain a good fit at cluster mass scales, we find 
that the Press-Schechter model overpredicts the number density of halos compared to 
the simulations by a weakly cosmology-dependent factor of 1.5-2 at galaxy and group 
masses. It is impossible to obtain a good fit over the entire range of masses simulated 
by adjusting 5c within reasonable bounds. 

Key words: large-scale structure of universe - dark matter - cosmology: theory - 
cosmic microwave background 



1 INTRODUCTION 

The COBE DMR det ection of anisotropi es in the cosmic mi- 
crowave background ( Smoot et al. 1992| ) made it very clear 
that the 'standard' structure formation scenario of c old dark 
matter (Blumenthal et al. 1984; Davis et al. 1981:) cannot 



simultaneously account for fluctuations on very large and 
very small scales. That model made several very restrictive 
assumptions about cosmological parameters - that space- 
time is homogenous, isotropic and globally flat; that there is 
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no cosmological constant; that fluctuations from homogene- 
ity are Gaussian-distributed and nearly scale-independent 
at horizon crossing; that the Hubble parameter h = Hq / 
(100 km s~^ Mpc~^) is 0.5; and that the number of free pa- 
rameters is minimized. The obvious fixes to the problem of 
excess small-scale power (when normalizing power spectra 
to the COBE anisotropy) are to make one of the following 
modifications to the model: 

(i) tilt the primordial spectrum, 

(ii) allow a nonzero cosmological constant but retain glob- 
ally fiat geometry, 

(iii) allow the universe to be open, 

(iv) add hot dark matter (i.e., neutrinos with masses of a 
few eV), or 



© 0000 RAS 



M. A. K. Gross et al. 



(v) lower the Hubble parameter much further {h ~ 0.3- 
0.4). 

Each of these modifications adds only one free parameter to 
the cosmology. In this paper, we consider the most viable 
models from each class above except the last, and simulate 
them with an A''-body code in two suites, with equivalent 
initial conditions across all the models. We do not consider a 



low-ii'o' model (Bartlett et al. 1995) because of increasingly 



solid observational evidence that h > 0.5. 

Deciding on cosmological parameters is to some extent 
an iterativ e pro cess. Much can be done using the Press- 
Schechter (1974) approximation, but the assumptions that 



go into it are n ot necessarily realistic (for exa mple, spherica l 
symmetry - see Jain & Bertschinger 1994 and Monaco 1995). 



Therefore, it is useful as a first approximation to calculating 
the mass functions, and we use it to perform an approxi- 
mate cluster normalization, using guesses about other cos- 
mological parameters. We run a set of simulations and use 
them to test the Press-Schechter approximation, and make 
several preliminary comparisons to ob servational data. In 
a companion paper (Gross et al. 1998), we recalibrate the 
Press-Schechter approximation and use it to derive refined 
estimates of model normalization and flo from several differ- 
ent data sets, and make more careful comparisons to cluster 
abundance. Subsequent papers will use simulations based on 
the refined normalizations. 

we describe our specific models from each 



In section 2.1 



class of CDM-variant models an d ex plain why we chose the 
parameters as we did. In section 2.2, we briefly describe the 



implementation of the particle-mesh algorithm we used for 
this study. We explain our halo finding algorithm and the 
effect of mass resolution upon it in section 3^ and report the 
simulation results in section Ml. Finally, in section B, we give 
our conclusions. 



2 SIMULATIONS 

2.1 Models 

Given the long list of modifications to the cold dark mat- 
ter scenario in the previous section, we could construct a 
model by adjusting every parameter in order to fit all the 
available observational data. However, in addition to being 
aesthetically displeasing, the physical significance of such a 
model would be unclear. As a result, we have tried to min- 
imize the number of modifications to the relatively simple 
standard cold dark matter scenario by investigating each of 
the modifications mentioned above in a separate model. The 
exceptions to this policy are that in addition to any one of 
modifications (ii)-(iv), we allow a small tilt, up to n = 0.9, 
in order to simultaneously fit the COBE and cluster data, 
and we allow the Hubble parameter to be adjusted within 
reasonable observational bounds according to the require- 
ments of the model. Larger tilts are not allowed because 
they tend to cause disagreements with high-multipole cos- 
mic microwave background data. 

We explore the large parameter space by running a 
large suite of linear calculations and comparing the output 
to appropriate observational constraints. Constraints that 
we consider in choosing model parameters for more detailed 
nonlinear analysis are: 



(i) the abundance of Abell clusters, as measured by X- 
ray temperature pr ofiles (White, Efstathiou, & Frenk 199J 
Biviano et al. 19931, hereafter WEF93 and BGGMM93, re- 



spectively). We assume that cluster masses may be underes- 
timated by up to a factor of two, motivated by results from 



cluster dens i ty mapping with gr a vitational lensing ( Squire! 
et al. 1996^ ^quires et al. 1997 J; |Miralda-Escude & Babul 



1995| ; |Wu fc Fang 1996| ; |Wu fc Fang 1997] ; figu: 



ire 



(ii) microwave background anisotropics for (. <. 800 (fig- 
ure 0) a s measured by s everal recent CMB detection exper- 



imen t s (Tegmark 1996| ; 



Netterfie ld et al. 
199^ ; [Piatt et al. 1997^ figure g) 



1997 



Scott et al 



(iii) 'bulk fiow' peculiar velocity meas urements and re- 
sulting constraints on the power spectrum ( Dekel et al. 1997 ; 
Kolatt fc Dekel 1997 ; figure ]^ . The linear estimates of these 
parameters are shown in table |l| and in figures |l|, y, and ^ 
for the models we consider. 

In most of the previous work with modified CDM mod- 
els, the most 'extreme' values of the model parameters have 
been chosen (i.e. as far from SCDM as was considered ob- 
servationally plausible). For example, low-fio models typi- 
cally have values of Qq ~ 0.2 — 0.3. However, in every case, 
while solving some of the problems with SCDM, this in- 
troduces new problems or confiicts with other observational 
constraints. Thus our approach will be somewhat different. 
We use our previous experience with linear and non-linear 
tests of CDM-variant models, as well as the published results 
of others, to find models that represent a 'middle ground' 
between SCDM and the most extreme version of the partic- 
ular class of model. In this way, we hope to choose the 'best' 
rather than the most extreme case, and to identify models 
that agree with the widest possible range of observations. 

For most models, we presume a baryon abundance of 
Sib = 0.025/i~^, consistent with the Tytler, Fan, & Buries 

measurement .n Nor- 



(199m) cosmic deuterium abundance 



malization is accomplished by calculating low mul tipoles 



using an e nhanced version of the linear code from Holtz- 
man (198!: ) and comparin g to the four-year COBE DMR 
anisotropy measurements ( Gorski et al. 1996 ; Gorski, pri- 
vate communication). 

For comparison to other studies, we also simulated the 
standard cold dark matter (SCDM) model with bias b = 
ag^ = 1.5. That model is intended to approximately match 
observed cluster abundances at the cost of being inconsistent 
with the COBE anisotropy measurements. For this model, 
we presumed there were no bary ons in the Universe , and 
used the BBKS transfer function (Bardeen et al. 1986) used 
in previous studies, that is. 



^ ' (2.34g)2 

(1 + 3.89q + {m.lqf + {5A6qf + (6.715)" 



-1/2 



(1) 



Buries & Tytler (1997a, 1997b) have very recently remeasured 



the deuterium abundance and found it to be 20 per cent lower, 
0.019 ± 0.001. This makes a very small change in the power 
spectrum, and the most significant effect is to make agreement 
with high-£ cosmic microwave background measurements (fig- 
ure U) more difficult. The height of the first Dopplcr peak depends 
strongly on Q\^. 
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with q — kh~^ and A adjusted so that the rms fractional 
variance in mass in spheres of radius 8 h~^ Mpc estimated 
using Unear theory is as = 0.667. 

The simplest way to solve the problem of excess power 
on small scales is by 'tilting' the spectrum, that is, by chang- 
ing the Ak factor in equation (hi) to ^fc", with n < 1. How- 
ever, the price paid is that choosing n < 1 also reduces the 
amplitude of the first 'Doppler peak' in the small-angle cos- 
mic microwave background spectrum. We find that n ~ 0.9 
is the largest allowable tilt that is still marginally consis- 
tent with the large-multipole cosmic microwave background 



data. "W 



ith n — 0.0. when wc COBE normaligc the model wc 



find th at it tends to overproduce clusters at Ad — 6 x fO^ 



h~^ Mq according to the Press-Schechter estimate, unless 
we use a rather low value of the Hubble parameter, h = 0.45. 
Although this is not favored by most of the current obser- 
vational data, we conclude that this choice of parameters 
consitutes the best compromise amongst the observational 
constraints that we have imposed. 

Another fix is to add a little hot dark matter, usually 
assumed to be in the form of a massive neutrino. Previ- 
ously studied versions of this class of model typically postu- 



late a tingle species of neutrino with significant mass, and 
a fraction iz^ — u.6 ot the critical density m the torm ot 
hot dark matter. This model was ruled out, based on its 
inability to reproduce the observed ab undances of Damped 
Lyman a systems (DLAS) at z ~ 3 (Kauffmann & Char- 



lot lOOJu Klypin et al. 1995a). Models with lower fractions 



of hot dark matter {Q^ = 0.2) a re more plausibly c onsis- 
tent with constraints from DLAS (Klypin et al. 1995a), but 
still have too much small scale power and thus overproduce 



cluste r^ at z = 0. However, as pointed out by Primack et al 



(1995; see also Pogosyan & Starobmsky lyyb ) if the hot dark 



matter is divided into two species of neutrino with equal 
masses, the power on cluster scales is reduced by 20 per 
cent without affecting smaller or larger scales. This lowers 
cluster abundances without worsening potential early struc- 
ture formation problems (small-scale power) or compromis- 
ing the COBE normalization. We find reasonable agreement 
with observed cluster abundances with Qi, = 0.2, A^^ — 2, 
h = 0.50, and n = 1; or alternatively, with a small tilt 
(n — 0.9) and a higher Hubble parameter {h — 0.6). We 
chose the former, based on concerns about the age of the 
Universe. However, we ran the simulations before the Hip- 



(ilo = 0.4) and a tilt (n = 0.9) to reduce small-scale power 
and correlations. 

Many observers favor an open cosmology and a high 
Hubble parameter, consistent with local density estimates 
and the Hubble Key Project. The lowest reasonable value 
of Jlo, given initial Gaussian fluctuations as assumed in 
all CDM-variant models considered here, is constrained to 



be above 0.3 at > 4g confidence (Nusser fc Dekel 199^ 



Dekel & Rees 1994 



Bernardeau et al. 1995). We 



cf. also 

adopt ilo ~ 0.5 as a 'reasonable' value for OCDM, noting 
that even this relatively high f2o leads to a power spectrum 
lower than that indicated by the potent analysis ( Kolatt & 



Dekel 1997; see also figure W). Our linear code is not capable 



of determining low multipole cosmic microwave background 
fluctuations for OCDM, as it uses a plane wave expansion 
that is only appropriate for flat cosmologies. Instead, we 
use fitting functions for the n ormalization ShJ^o ) and the 
transfer function T{k) given by Liddle et al. (199q , hereafter 
LLRV96).| 

Figure bJ summarizes the expected mass functions on 
the group and cluster mass scales, as estimated from the 
Press-Schechter approximation, with a Gaussian fi lter. U s- 



ing the calibr ation with A-body simulations from Borgani 
(1997q), we use Sc.s = 1.5 for the model with mas- 



ct al. 



sive neutrinos and 5c, g — 1.3 for all other models, in this 
figure (but cf. Table W for best-fit (5c, g and 5c,t to our simula- 
tion results) . The observational cluster abundance estimates 
plotted are in reasonable agreement with these mass func- 
tions, especially if the mass estimates are low as indicated 
by some gravitational lensing estimates. 

In figure H, we compare each model to several recent 
CMB measurements , using the CMBFAST program of Seljak 



Zaldarriaga (1996). We also show the four most recently 



parc os recalibration of the age of globular clu sters (Reid 
19971; Idratton et al. 19971; IChaboyer et al. 19971). This con- 



straint has now been considerably weakened, and h — 0.6 
or even 0.65 would lead to ages consistent with the present 
estimates. 



announced CMB results on the figure. Not shown are sys- 
tematic calibration errors of 14 per cent for Saskatoon and 
20 per cent for Python HI. Note that the OCDM model 
is strongly inconsistent with the Saskatoon points, and our 
choice of Jlo = 0.5 is at the 95 per cent confidence lower limit 
for an open model (Lineweaver & Barbosa 1997). Also note 
that the models with even the relatively mild tilt of n = 0.9 
are at best in marginal agreement with the Saskatoon data 
around the first Doppler peak. 

Figure hi shows the linear power spectra at the present 
epoch. As one might expect, all the spectra nearly cross at 
a wavenumber of a few tenths h Mpc~^, corresponding to 
cluster scales. Also, we show some of the window functions 
used in the normalization procedure described above. 



Th3 Qa 7^ class of models has been well studied, typ- 



After we r an thi s model. LLRV96 w as superseded by 



ically with S2o = 0.3 and h ~ 0.65-0.7. However, analysis 
of earlier A'^-body simulations has shown that when non- 
linear effects are included, this model produces a power spec- 
trum/correlation function with too high an amplitude on 
small spatial scales compared to observations, unless galax- 



White (1997) and Hu & White (1997). Those papers' erg values 



Bunn & 



les arc rstr 



iliuiigl.y diiLi-bidSL'd with luspuLt Lu Llie daik iiiaL- 



agrec to high precision wit h LLRV96 if one lowers f^b from 0.025 
h-^ to 0.015 h'2, which ^unn fc White (1997 ) favor anyway. 
However, the transfer function shapes are somewhat difEerent, 
and the LLRV96 normalization is to the COBE 2-year data, so 
the power on scales of a few hundred fc~^ kpc may be up to 20 
per ce nt low compared to Bunn & White (1997) and Hu & White 



ter (Kl ypin, Primac k ._ fc Holtzman 1006 , hereafter KPH06 



(1997). Using a BBKS-style fit as all three papers do, rather than 



Jenkins et al. 1997 ). phigna et al. (1997 ) have also shown 
that the void probability function for this model is in dis- 
agreement with observations. Therefore we have chosen a 
model with a slightly higher value of the matter density 



integrating the Boltzmann equation directly, introduces an error 
of similar m agnitude, even with t he improved shape parameter 
described in Hu & Sugiyama (199f:, equation D-29). We the refore 
neglect the diff'erence between the Bunn Sz White (1997) and 
LLRV96 spectra. 
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Table 1. Model parameters and linear results for both simulation suites. 



Model 


Age'' 


h" 


Oo 


f^c 


Qb 


n^ 


f^A 


n" 


N^ 


^s 


~4 


''so 


^5 


observations 






















375 


5 X 10-6 


1-cr errors 
























85 


2 X 10-'' 


CHDM-2i/ 


13.0 


0.5 


1.0 


0.7 


0.1 


0.2 


0.0 


1.0 


2 


0.719 


0.719 


399 


6 X lO-fS 


OCDM 


12.3 


0.6 


0.5 


0.431 


0.069 


0. 


0.0 


1.0 





0.773 


0.581 


254 


3 X lO^fS 


SCDM 


13.0 


0.5 


1.0 


1.0 


0.0 


0.0 


0.0 


1.0 





0.667 


0.667 


192 


2 X lO-f* 


TCDM 


14.5 


0.45 


1.0 


0.9 


0.1 


0.0 


0.0 


0.9 





0.732 


0.732 


270 


5 X lO-is 


TACDM 


14.5 


0.6 


0.4 


0.365 


0.035 


0.0 


0.6 


0.9 





0.878 


0.572 


335 


2 X 10-6 



Mpc" 



n^h-' 



Time since the Big Bang in Gyr. 

Presumed Hubble parameter, in units of 100 km s 

'Tilt' of the primordial spectrum; P{k) oc k". 

Number of massive neutrinos presumed. The equivalent mass of a neutrino is m^ 

rms mass fluctuation in a sphere of radius 8 h~^ Mpc 
~ — oO-46-O.lOno 

(1996) calculate erg = 0.52 ± 0.04 in order to fit cluster temperatures assuming /3 = (KE)^ 
However, new simulations (Frenk ot al. 1998) show that /3 sa 1.17. co rresponding to ag 



■ 92 eV. 



f ^ 



c!g,QQ ' ° for A = models, and rTg,Q,Q ' ° for SIa + f^O = 1 models. |Eke, Cole, 



Frenk 



0.61 ±0.05. 
s rms velocity in a sphere of radius 50 h~^ Mpc (Dekel et al. 1997). Note that the observational value is 
for one particular 50 h~^ Mpc sphere around the Local Group, and the simulation values are the rms value 
for a distribution of randomly placed spheres. These are not the same, so we can't use the observations to 
rigorously define confidence limits for the simulation quantities. 

'' Estimated number density of clusters of mass > 6 X 10^* h~^ /V/"g. in h'^ Mnc-^. from Press-Sr.hechter 
theory with a Gaussian window functi on. S,- a is 1.5 for G HDM models ( Walter fc Klypin 1996 ; Borgani et al. 
19971) and 1.3 for all other models (Liddle et al. 1996, KPH96). Masses near the center of the allowable 



range for cluster data (WEF93, BGGMM93) are used. Note that uncertainties in the masses of measured 
clusters mean that the masses for which densities were measured could shift coherently up to a factor of two 
above the reported value of 4.2 X lO^** h~^ A^q. This is a naive app roach which we us e only for our first 
iteration of model parameters. Wc use more sophisticated methods in Glross et al. 



In figure kl we compare our models to the mat ter power 
spectru m re cently measured from bulk flows by Kolatt & 
Dekel (1997). We only use the three data points that Ko- 
latt & Dekel use for their own statistical analysis, because 
for larger wavenumbers smoothing lowers the power signif- 
icantly. SCDM disagrees at about the 2.5a" level, also re- 
flected in its low value of V50 in table |I[ OCDM and TACDM 
disagree because the value of P(~ 0.1 h Mpc-^) is fixed by 
comparin g the observed densi ty of clusters (WEF93, BG 



GMM93, [Borgani et al. 1997aD to the Press-Schechter pre- 
diction, and they have low values of f{flo, SIa) = Da/ Da ~ 
f2o'6^ where D{Q,f,,K, t) is the linear growth factor and a{t) is 
the expansion parameter. The combination of cluster abun- 
dances and bulk-fiow power spectrum measurements favors 
/ ~ 1, for the currently favored classes of CDM- variant mod- 
els. 

There is currently significant controversy over the 
proper normalization of models, and our OCDM and 
TACDM normalizations are higher than the recent fits re- 
ported in Eke, Cole, fc Frenk (199q), based on clu ster X-ray 
temperature distributions (Henry & Arnaud 1991 ),n though 



they are consistent with the older analyses of WEF93 and 



Bor 



the ne wer cluster velocity dis persion measuremei its of 

gani et a.1. (1997a ). Pen (1996) has reanalyzed the 

& Fren : (1996) calculation, and he gets slightly higher low- 



Eke, Cole 



Note that 



Eke Cole, fc Frenk fl 99f ) discovered two compensat- 
thc Henry & Arnaud (1991) analysis; an arithmetical 
error of a factor of about 4 



ing errors in 

error of a factor of 4.2 and a binnin, 

in the o ther dir ection. The errors have also been noted in |Viana 

& Liddl ; (1996| ) 



Qo normalizations of erg = 0.86 and 0.72 for our TACDM 
and OCDM models, respectively. These normalizations are 
close to those we have chosen (table |l|) . 



2.2 Algorithm 

A classic problem with gravitational simulations is the 'over- 
merging' problem, where small scale structure in highly over- 
dense regions is not resolved. Part of the problem is physical 
- real galaxies form much denser cores than dissipationless 
hal os can, because the baryons can dis sipate energy (but 
cf. Klypin, Gottlober, & Kravtsov 1997). Aside from that. 



numerical limitations can make the problem vastly worse. 
There are two numerical effects to consider: force resolution 
and sampling of initial conditions and bound structures. Im- 
proving either of these requires vast amounts of memory and 
processing time, so there is an inherent tradeoff. 

Recently, the more popular approach has been to im- 



prove the forces by using hy brid (Hockney 

Couchman 199l|; |Xu 199. 

(|Kravtsov, Klypin, & Khokhlov 1997 



: Eastwood 198J 



i| 

q, for exampl e) or adaptive-mesh 

for example) force 



solvers, at the expense of either poor sampling of initial 
fiuctuations or small box sizes. We choose a complementary 
approach, where we try to balance the sampling of density in 
a large box with the force resolution. We still require a large 
dynamic range in order to sample small scales well and si- 
multaneously simulate a large volume for comparison to red- 
shift surveys. Since the two requirements imply an enormous 
number of particles, computer time limitations force us to 
use th e fastest code available. We ch oose a standard particle- 
mesh (Hockney fc Eastwood 1988) algorithm, parallelized 
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Figure 1. Expected mass functions for all models estimated from 
the Press-Schechter approximation with a Gaussian filter, using 
5c, g = 1-5 for the model with massive neutrinos and 5c, g = 1-3 
for all other models. The two data points shown correspond 
to observational estimates of cluster abundance (WEF93, BG- 
GMM93 1. Note that cluste r densitv mapping via gravitational 



GMM93 1. Note that cluste r densitv mapping ' 
lensing (pauircs ct al. 199fi|: -jauircs ct al. 1997 ; 
& Babu 1995 ' Wu & Fang 1997 ) may indicate tl 



Miralda-Escude 



Wu & Fang 1997 ) may indicate that X-ray masses 



are systematically low, and the masses can plausibly be raised by 
a factor of up to two, which corresponds to the horizontal line on 
the right of each cluster data point. 



to run on a distributed-memory message-passing system.^ 
This type of code produces adequate forces at about 1.5 grid 
cells (KNP97, Appendix A), but we double this distance to 
be conservative. So, we require that we have 3'^ times as 
many grid cells as particles, for the high resolution suite. 
We choose a grid cell size of 65 h~^ kpc, with A'^g — 1152'^ 
grid cells and Np = 384'^ = 57 million 'cold' particles. For 
the large-volume case, we wish to follow the dynamics only 
of clusters of galaxies, so we can afford to coarsen the den- 
sity grid slightly. We find that a cell size of 390 h~^ kpc is 
adequate for following the dynamics of J> 10^^ h~^ Mq ob- 
jects, and expect information about smaller objects to come 
from the high resolution simulations. The slight coarsening 
of the density resulted in a substantial advantage in running 
time. 

Initial condi tions were calculated using a parallelized 



Zel'dovich (197C) approximation. For CHDM models, we 



started with a uniform grid of cold particles, and two neu- 
trinos at the position of every cold particle. Cold parti- 
cles and neutrinos were offset from the grid using separate 
cold-l-baryon and hot power spectra, and consistent veloci- 
ties were derived from the offsets using scale-dependent lin- 



* Specifically, the Cornell Theory Center SP2, but the code is 
portable to any system supporting MPI, including heterogeneous 
workstation clusters and most modern supercomputers. 



ear gr owth rates calculated by a refinement of the Holtzman 
(19891) code. In addition, equal and opposite random ther- 



mal velocities were chosen for each pair of neutrinos from a 



redshifted relativistic Fermi-Dirac distribution ( Klypin et al 
1993| ). 

We adopted the form of the equations of motion used 



Kates, Kotok, & Klypin (1991) generalized to arbitrary 



cosmology: 



v= 




3 

~ 2a 


5p 

f^Pc' 


dp 

da 


- = — Q 


V<^(xi), 




da 


Pi 




he 


Friedmann 


equation wit 



(2) 



(3) 



(4) 



able Hat, 

a — a~ ' V ^c + 0,^ -\- Q,t^a^ + Q,]^a. (5) 

Time discretization was a standard 'leapfrog' scheme (cf. 



Hockney & Eastwood 198S), with even steps in the expan- 
sion parameter a. To reduce the expense of the simulations, 
the timestep was chosen only to stabilize bound structures 
at the final timestep, rather than keep all structures on 
the scale of the grid spacing stable. This is only a problem 
for the cores of clusters, which have the highest velocities. 
For clusters, we presume an upper bound of particle veloci- 
ties of 1200 km s~^ today and a minimum diameter of any 
given bound structure equal to the linear cell size. Stability 
for such an object requires that particles take at least one 
timestep to traverse the object. So, the required condition 



Aa = dAi < 



HoL 



A. 



1/3 



(6) 



independent of cosmology because the condition is evaluated 
at the present epoch and HoL is chosen to be the same for 
all models. Plugging in Umax = 1200 km s~^, HqL — 7500 
km s"^ and A^g = 1152^ gives Aa ^ 0.005, or 200 timesteps 
for the high resolution suite. Such a low Umax will not model 
the interiors of clusters well, since they are observed to have 
velocity dispersions larger than that, but to remain bound 
to the cluster, particles have the much looser requirement 
that they not traverse the whole cluster in one timestep. As 
large cluster radii are up to about 50 grid cells, the effective 
stability limit is 20 per cent the speed of light inside a large 
cluster, for the high resolution suite, presuming that the 
cluster is adequately modeled by an isothermal sphere. We 
checked that the choice of timestep was adequate by running 
a 25 /i~^ Mpc box CHDM-2j^ simulation with 384^ grid cells 
(which has the same 65 h~^ kpc cell size as the high reso- 
lution suite) for 200 timesteps and for 300 timesteps. The 
resulting mass functions were not significantly different. For 
the large volume suite, clusters do not cover nearly as many 
cells as in the high resolution suite, and so the velocity limit 
is much higher. Our choice of 150 timesteps corresponds to 
a limiting speed of 5000 km s~^ if particles are not to cross 
one cell in a timestep. The suite parameters are summarized 
in table |.| 

^ Though the implementation, and especially parallelization, of 
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Comparison of cosmological models to selected CMB observations 



100 



a. 




Figure 2. Model C omparison to Cosmic Microwave Background. All models except SCDM are consistent with COBE four year data 
( Gorski ct al. 199i: : Gors ki. private communicatio n). Circles, solid squa res, open squares and asterisks ar e the COBE four year power 
spectrum ( Tegmark 1996| ), Saskatoon 1995 results ( |Nctterfield et al. 1997 ), CAT detection ( ^cott et al. 1996| ) and Python III results (Piatt 
et al. I9M), respectively. Not shown are systematic normaliz ation errors of 14 and 20 per cent, for Saskatoon and Python III, respectively. 
The curves are all calculated using the CMBFAST program of peljak fc Zaldarriaga (1996 ). Cosmological parameters correspond to models 
considered in this paper, except for SCDM. The normalization is adjusted so that the low harmonics match the output of our linear 
code. CMBFAST is capable of calculating larger multipoles than our linear code. SCDM is shown here with Oj, = 0.1, since all the high-£ 
features in the CMB spectrum are dependent upon baryon interactions, but was actually simulated with no baryons. 



The Zel'dovicli (197C) approximation is only valid when 



the rms fluctuations are much less than 1. In practice, one 
picks a starting time early enough so that linear theory 



the two-species particle mesh code described above is much less 
trivial than one might suppose, discussion of the code has been 
omitted for space considerations. The interested reader may find a 



detailed description of the code , its impleme ntation on the Cornell 
SP2, and several code tests in pross (1997]). 
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Table 2. Simulation parameters for both simulation suites 
Suite Box size N^^n^ 



Box size 
h^^ Mpc 



Cell size 
h^^ kpc 



cold 



N. 



(He + ^b) h-^ Mq 



steps 



mm •' 



Co h-^ M(r 



high res 


75 


11523 


65 


384^ 


2.09 X 10^ 


200 


3.4 X 10" 


0.078 


low res 


300 


768^ 


390 


384^ 


1.34 X 10" 


150 


7.3 X lO^^ 


0.043 



" Number of cold particles; for models with massive neutrinos, A'l,^^ = 2N^^i^. 
^ Mass of cold particles; for models with massive neutrinos, Mjiot = A/(.oidf!i//2(f2c + ^h)- 
'^ Halo detection cutoff, from the restriction that halos must be larger than the grid size. 
'' Fractional error in mass for the smallest halos identifiable. 



Linear power spectra 



Linear power spectrmn comparision to l^iilk flow nieasnrements 




10-^ 10 



10-2 10-' 10° 

A: (h Mpc-i) 



CHDM-2;^ - 

OCDM - 

SCDMi)= 1.5 - 

TCDM - 

TACDM - 

Kolatt et al. 1997 f 



5^ 




A- (h Mpc-') 



Figure 3. Linear power spectra used in our simulation suites. 
Also shown are two of the window functions used in nor- 
malizing the models: k^W'^{rk) with r = 8 h~^ Mpc for 
cr| and W^{rk) with r = 50 h'^ Mpc for V^g. Here 
W{x) = 3 [sin(x) — xcos(a;)] a::~3. Also shown (for illustrative 
purposes only) is the equivalent window function for approxi- 
mate COBE normalization using the pure Sachs-Wolfe effect, 
Jlo('^hfc)/(27rd^fc'^) where jio{x) is the 10th order spherical Bessel 
function and dji is the horizon distance. The version plotted 
has the amplitude raised by a factor of 100 for visibility and 
uses dji = 2c/Hq = 6000 h~^ Mpc, which is appropriate for 
Qo = 1- For OCDM, the horizon distance is 7470 h~^ Mpc and 
for TACDM, it is 8810 h~^ Mpc, so the window function moves 
a small distance to smaller k in those cases. A similar window 
function for cluster abundance doesn't exist because it doesn't 
have the form of a convolution. In an extremely rough sense, the 
scales are comparable to those sampled by eg. 



brings the rms fluctuations weU below 1. The initial time 
was chosen so that the rms overdensity on the grid scale was 
5rms <j 0.2. This was z = 30-60, depending on the model. 
Particle data and halo catalogs were stored at four equally 
spaced intervals in a during execution. The large volume 
simulation suite used the same starting times as the high 
resolution suite even though they could have been started 
somewhat later due to the poorer resolution. The extra com- 



Figure 4. Linear power spectrum comparison to bulk-flow mea- 
surements. The curves are all a magnification of figure M, multi- 
plicd bv f^ fOn . n A ) = (aD/aD)^. The three data points are from 
Kolatt &: Dckcl (1997]). f(^n.^\ ) was calculated exactly, using 
equation (C.3.14) of Gross (1997) and its analytic derivative. 



putation involved is about one timestep and is therefore neg- 
ligible. 

Random numbers are necessary to model inflation- 
generated Gaussian fluctuations and random phases in the 
density fleld. Such randomness introduces highly significant 
variation from simulation to simulation, commonly referred 
to as 'cosmic variance.' Because we can only observe one 
universe, quantifying the effect of cosmic variance is very 
important and is a separate issue from variations between 
models due to different physics. In these suites, we have sep- 
arated the effects by picking a single random number seed 
for each suite, checking that the largest 26 waves do not have 
any fluctuations larger than a factor of 2, and rerunning one 
model with a different seed. That is, within each suite, the 
random numbers for each model within a suite are all the 
same, and large wavelength fluctuations are restricted to a 
smaller range than Gaussian statistics would permit, in an 
attempt to prevent rare statistical fluke s from compromisin g 
expensive simulations (as happened in |Klypin et al. 1993 ). 
This means the structures are approximately in the same 
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place, 

criterion discussed in section |2.1 



and when one also considers the cluster abundance 
there is roughly the same 



number, distribution and positions of 5 x 10 h~ Mq clus- 
ters in all the models in a given suite. Note that the models 
do have different power spectra and fluctuation growth rates, 
so distributions can differ for objects with different masses. 



3 DARK MATTER HALOS 

3.1 Halo finding algorithm 

We identify dark matter halos using a spherical overden- 
sity algorithm similar to that of KNP97, with some of the 
limitations removed: 

(i) We define candidate halos as the centers of all den- 
sity maxima containing an overdensity greater than S = 
5p/Qopc = 50. A density maximum is defined as a cell whose 
density is greater than its six Cartesian neighbors. Just in 
case there are other halos hiding in those six neighbors, we 
also consider each of them to be candidates. Note that the 
finite grid size (as in all other grid-ba sed halo finders, such as 
DENMAX, Glelb & Bertschinger 1994) will introduce a mini- 



mum separation between halos, which may cause small halos 
to be missed, which in turn will require a mass cut. 

(ii) Each candidate halo then has the location of its center 
set iteratively to the center of mass of all the particles inside 
a sphere of diameter equal to the cell size (65 h~^ kpc in our 
high resolution case). Halos are expected to have a minimum 
size of the order of the grid size, so this procedure moves 
the candidate halo to the peak of the density maximum. Of 
course, since we have defined more than one candidate for 
each detected maximum, some candidates will converge on 
the same halo. The smaller mass object in a given pair is 
removed if the distance between the centers of mass is less 
than half the grid spacing. 

(iii) We perform a central overdensity cut. All halos that 
don't enclose a mean overdensity sufficient for virialization 
according to the spherical collapse model (see Gross 1997 
appendix C, and references therein) at the end of the center- 
of-mass detection phase are presumed not to be virialized 
objects and are discarded. Typically, this reduces the num- 
ber of halos by a factor of 2-3, though the number is model 
dependent. 

(iv) We now estimate at what radius the mean enclosed 
overdensity S = Sp/ilpc falls to S^ii, the virial radius of the 
halo in spherical infall models.n For each halo, we count the 
number of particles within five radii up to five grid cells (325 
h~^ kpc in this case) away from the center, convert that to 
density, and interpolate the radius at which S — S^ir (^vir) 
using power-law cubic splines. If five radii is not large enough 
to enclose rvir, we search five more radii, each twice as long 
as the original radii. This is repeated until we enclose rvir. 

(v) We define the mass of the halo as the mass enclosed 
in Tvir. The velocity is the mean velocity of all the particles 
within rvir. 



(vi) In general, the largest halos in a high resolution run 
contain much resolved but bound substructure. Because we 
search for the S — S^h radius, we detect the same regions 
of space dozens of times for the largest halos. To remove 
'double-counted' halos, the halos are searched in reverse or- 
der by mass to see if they enclose the centers of any smaller 
halos. If so, the smaller halo is thrown away. Note that 
the ordering is important because three-body intersections 
would be non-deterministic otherwise, and throwing away 
halos that only intersect is too stringent. 

One limitation of this algorithm is that it presumes all 
halos are spherically symmetric, which is demonstrably un- 
true. But the effect on the mass function is random, rather 
than systematic, and finding the halos with an algorithm 
generalized to ellipsoidal distributions does not change the 
mass function significantly, even though it changes the pa- 
rameters of individual halos. Because the halos have finite 
size, one cannot perform mass-weighted correlation func- 
tion analyses, for distances less than the largest halo radius 
(about 2-3 h~^ Mpc in radius, typically). 

The other limitation is the use of the density grid to 
identify halo candidates. If one considers a worst-case iden- 
tification where a large number of particles all collect in one 
corner of a grid cell, in order to guarantee that all nearby 
halos are identified, one must draw a sphere which encloses 
the entire cell, of radius 



= ^L/iVg'/^ 



(7) 



^ Not e that our definition of iSvir is related to Eke, Cole, & Frenk 
(1996) by 1 + <5vir = Aecf/Q. Our choice is appropriate for the 



density field calculations in an A'^-body code. 



where L is the length of one side of the computational vol- 
ume and TVg is the number of grid cells. If halos happen to 
be bigger than that, then the last step of the halo catalog 
generator makes it unimportant that we couldn't see nearby 
structure. Fortunately, halo extent is trivially related to halo 
mass because we have defined both where the mean overden- 
sity is 5 = 5vir. 



3.2 The effect of mass resolution 

To what extent should you, the reader, trust the mass func- 
tions presented in this paper? To answer that, one must con- 
sider several effects. A typical feature in a mass function is 
that the large-mass end becomes 'wiggly,' usually blamed on 
the scarcity of high mass halos combined with cosmic vari- 
ance. There is a related effect at somewhat smaller masses, 
since very large halos tend to have somewhat massive com- 
panions. For example, in most models in our high resolution 
suite, 5 X 10^'' h~^ Mq objects are fairly rare, but it is com- 
mon to see them as companions for 10^^ h~^ Mq objects. 
So, the wiggles may propagate down the mass function, and 
cosmic variance may have a significant effect on more than 
just the highest mass scales. 

Cosmic variance fortunately leaves a signature, in that 
the mass function is not smooth at high masses. But, it 
is quite important to figure out the limiting factors at low 
mass, where typical mass functions are quite smooth. What 
limits accuracy here are the effects of finite sized grids and 
finite numbers of particles. 

The effect of the finite sized grid in identifying maxima 
in the final particle distribution was discussed above, and 
one must merely translate the minimum radius of a halo 
rmin to a minimum mass. Since the halo radius and mass 
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are defined as enclosing a mean overdensity of Svir , the mass 
Mvir of a halo of radius r^ir is 



47t 3 

Mvir = (1 + c5vir)— f^OPcr-vir- 

So, a very conservative mass cut is 

M^in = (1 + ^vir)Y»0/Oc ^ ^ ^ ■ 



(8) 



(9) 



Plugging in values for the high resolution suite, the mass cut 
is 3.4 X lO^^rJo h~^ Mq. For simplicity, we make the same 
mass cut on all models, corresponding to fio = 1. 

One might worry that the central density cut described 
in the previous section could cut too many small halos, be- 
cause the fairly long timesteps used cause the density within 
the 'half-mass' radius to go down by ab out a factor of tw o 
if the timestep equals the stability limit (Quinn et al. 1997). 
We perform the central overdensity cut at r = L/2N^, 
but the proximity restriction used in deriving equation (bl) 
requires that halo radii in the final catalogs be at least 
\/3 L/N^ . If halo profiles fall at leas t as fast as r~^ (whereas 
the Navarro, Frenk, & White 199i profile says it should be 



much steeper than that near the virial radius), then the den- 
sity fed into the central overdensity cut should be at least 
a factor of 2^/3 ~ 3.4 greater than the virial density. This 
more than offsets the density smoothing due to timestepping 
at the stability limit, so we neglect the effect of time steps 
in our mass resolution analysis. Note that our timesteps are 
only near the stability limit for virial radii near the detection 
limit - otherwise, a particle takes many timesteps to cross a 
halo. Therefore, lowered densities due to long timesteps are 
only a concern for the smallest detectable halos. 

One might also worry that the quality of the force law 
at scales approaching the grid scale would also result in re- 
duced central density. At the 1.7 grid cells proximity cutoff, 
the point-mass potential in our simulations is about 90 per 
cent of the correct GM/r value. With such a force law, the 
virial theorem requires that the density be also 10 per cent 
low, to maintain the same velocity dispersion. Thus, some of 
the smallest halos around the mass cut will not make it into 
the catalog. In practice, the density profiles for the smallest 
halos are considerably noisier than 10 per cent due to as- 
phericity and background particles, so we neglect the effect 
of an oversoftened force law. 

Particle discreteness may also affect the halo mass func- 
tion, because random fluctuations may affect the detection 
of some of the smallest halos. We consider here how much 
significance we need to make the expected number of halos 
missed fewer than the number of halos in the simulation. 
Because halo boundaries are defined where the mean over- 
density 5 is (5vir and the mass of a particle is 



Afp = f2oPc 



N^ 



(10) 



where A'^p is the number of cold particles in the simulation,n 



^ For simplicity, we consider the significance of models with only 
one particle mass. For CHDM-2i', a given mass will always be 
represented by more particles than in SCDM. So, one can get a 
conservative estimate of the significance by only considering the 
cold particles. 



the number of particles inside a halo of mass Mvir in a sim- 
ulation box of size L is 

For counting N particles within r, the random variation in 
number is a = \/N . Let us suppose there are Ah halos 
above a given mass, and we wish to detect them all. We 
presume that counting halos is a Gaussian process and state 
that the n-sigma uncertainty in the detection of the halos 
corresponds to incorrectly detecting or missing a fraction 
erfc(n/\/2) of the halos. We require detection of all halos, 
so the fraction missed should be less than 1/ p^^L , where p is 
the number density of halos above the mass cutoff, and L^ is 
the volume of the simulation box where halos are identified. 
Inverting, we need detections of \/2erfc~'^ (l/phL"^) sigma. 
The density ph should really come from the Press-Schechter 
approximation, given a desired mass cutoff, but the inverse 
complementary error function is extremely insensitive to the 
value of its argument, once it becomes much less than one. 
As an example, for ph = 1 h'^ Mpc~"^ (appropriate for a mass 
cutoff a little below lO'^^ h~^ Mq for most models), we need 
to have at least 4.7a detections of all halos. Less significant 
detections mean it is likely some of them have been missed 
by random fiuctuations. This means every halo must contain 
at least 23 particles. More generally. 



A, 



erfc 



1 



PhL-^ 



(12) 



for the rather liberal restriction that we only require detec- 
tion of the halo. 

As an alternative cutoff criterion, requiring a 10 per cent 
or less 1-a error in mass is a more stringent requirement, 
and every halo must have at least 100 particles, since mass 
is determined by counting particles within several radii. If 
one requires a fractional error of / for a minimum halo mass 



of A/„^ 



one needs at least 



N^ = 



floPcL 



(13) 



particles in the simulation. The parameters used, and the 
effective / they allow, are shown in table ti. Figure H shows 
that, with grid sizes and mean interparticle spacings of the 
order of those used in our suite, the effect of lowering either 
the grid size or the mean interparticle spacing by a factor of 
two does not significantly affect the mass function. For this 
test, we raised the threshold for halo candidate identification 
from (5 = 50 in one cell to (5 = 70 because one isolated cold 
particle in the high Ag case gives 1 -|- J = 51.2. 

To explicitly test the effect of grid sizes on our mass 
functions, we ran five small simulations of the CHDM-2!/ 
model with Ag = 192^ grid cells and Ap = 3 x 64^ parti- 
cles, with various-sized boxes. Though these simulations are 
too small to generate meaningful mass functions on their 
own, collectively their upper envelope does match the Press- 
Schechter formula reasonably well, for 5c, g — 1.2 with a 
Gaussian filter. Figure shows the five different mass func- 
tions. Also shown are lower mass cuts, determined for ev- 
ery model using equation (S). Above the mass cutoffs, every 
mass function agrees with the one for the next smaller box. 
Well below, the mass function slopes are not steep enough, 
but they agree with the neighboring curves for significant 
distances below the mass cuts, so it may be reasonable to 
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Figure 5. Effect of raising the number of particles or the number 
of grid cells by a factor of 8 in a very small CHDM-2i/ 10 h^^ Mpc 
simulation with Ng = 128^ and Np = 3 X 64^. The mass functions 
are not significantly different. A somewhat low mass cutoff of 10^^ 
h~^ Mq has been applied. The high resolution suite has a linear 
cell size that is slightly smaller than the A^g = 128'^ runs shown 
in this figure. 



extrapolate the mass function further. Every halo detected 
by the halo finder is represented in the figure, and the loca- 
tions of the lower mass cuts are indicated by vertical lines. 
This test could conceivably overproduce clusters because of 
the extremely poor force and mass resolutions in the largest 
volume run - a cell width is about the size of an Abell ra- 
dius. This Sc result does persist for much larger simulations, 
as di.sciissed below. 



4 RESULTS 

The connection between simulations and observations is still 
fairly uncertain, and the least well determined portion of it 
is the galaxy identification procedure. It is therefore help- 
ful to do as much analysis as one can using quantities that 
are insensitive to the deta ils of galaxy formation . Currently, 
only bulk flow motions ([Kolatt fc Dekel 1997) provide a 



meaningful matter power spectrum, but the large smooth- 
ing required means that the comparison is best made to the 
linear power spectrum (see figure M). When investigating 
quantities derived from observations of galaxies (as the vast 
majority of astronomical observations are) one is forced to 
make assumptions based on expections about the nature of 
galaxy bias, for example the usual expectation that galaxies 
are more clustered than the dark matter. Figure m shows 
nonlinear real-space dark matter power spectra for all our 
models , compared to the APM rea l-space galaxy power spec- 
trum ([Baugh fc Efstathiou 1994- The OCDM model re- 



Figure 6. Effect of grid size on mass functions. The curves repre- 
sent very small simulations of various sizes. The Press-Schechter 
model was tuned to match the cluster-scale part of the mass func- 
tion in the largest box. Note that the envelope mirrors the Press- 
Schechter curve reasonably well, but each individual mass func- 
tion has a power-law index that is too shallow. Mass functions are 
limited at the large-mass end by statistics - one simply runs out 
of enough space to create objects in — and on the small-mass end 
by some fraction of the halos becoming as small as two grid cells, 
which means it is not guaranteed that the halo can be resolved 
from its neighbors, particularly if they are also small halos. The 
vertical lines represent the lower limit in mass for each run, above 
which all halos can be detected. 



and it is very diffic ult to explain physically, especially on 



such large scales (c f. Yepes et al. 1997 



Kauffmann, Nussor 



Steinmetz 1997). Additional arguments against strongly 



quires significant antibiasing and the ACDM model requires 
even more. There is no evidence for such strong antibiasing. 



scale-dependent antibiasing are given in KPH96. 

The process of galaxy formation is not well understood, 
so one could argue that perhaps there is some mechanism 
that would give us strong antibiasing. We have created an 
extreme model for galaxy formation designed to produce 
as much antibias as possible (cf. KPH96). Everywhere in 
the density grid, if there is more than 2.1 x 10^ h~^ Mq 
in a grid cell, we presume one galaxy forms there. That 
mass corresponds to slightly more than the mass due to one 
isolated particle in the high-resolution SCDM and TCDM 
simulations (which have the most massive particles in the 
suite). Such a limit is necessary to prevent placing excess 
power in the voids due to vestiges of the initial grid there. 
This is a highly unreasonable model for galaxy formation, as 
it says that the density of >2 x 10^^ h~^ Mq galaxies in the 
core of the Coma cluster should be the same as in the local 
group, and this is clearly ruled out observationally. However, 
even though there is significant antibias on small scales, it 
is only visible at scales smaller than about k — 1 h Mpc~^ 
(see figure 8), whereas antibiasing is needed on scales larger 
than that in order for OCDM or TACDM to be consistent 
with the APM power spectrum. Note that a possible way 
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Figure 7. Nonlinear real-space dark matter power spe ctra, com- 
pared t( the AP M real-space galaxy power spectrum (Baugh & 



Efstathiiu 1994) of galaxy number-count fluctuations. The simu- 
lation p )wer spectra shown here are a composite of the high and 



low res olution auitca, where data from a niodcl'a high reaolution 
run is used at large k and low resolution data is used at small k. 
Two different high resolution runs of the SCDM case are shown as 
a guide to how large cosmic variance is. The power in the second 
SCDM realization is 20—30 per cent lower than that in the first 
realization for 0.3 < k < 1 h Mpc~^. The APM data are presum- 
ably biased with respect to the matter power spectrum, and yet 



the OC pM, TCDM, and TACDM cases require the APM data 
to be si gnificantly antibiaacd with rcapeet to the dark matter, 
with b2 ^ 0.6 for OCDM and TCDM and fe^ ^ o.5 for TACDM 
at fc ~ 1 h Mpc~^. If APM misses galaxies in clustered regions, 
that would give a low power spectrum on scales of fc > 1 /i Mpc~^ 
(see text). 



out of the antibiasing requirement is to note that the APM 
survey is incomplete in clustered regions, which will raise 
the 'true' power spectrum above the APM measurement on 
small scales (ZabludofT, private communication). 

The APM power spectrum is not the only power spec- 
trum that has been measured. However, to compare to other 
measurements, it is usually required to calculate model red- 
shift space power spectra. Going to redshift space signif- 
icantly reduces power on scales of interest, since typical 
dispersion velocities of 1000-2000 km s~^ in clusters cor- 
resp on d to a scatter of 10-20 h~^ Mpc in distance. In per- 



forminj this operation on particles, the power should be 
viewed as a lower bound, because there may be significant 
velocity bias (Carlberg, Couchman, fc Thomas 199C; Sum 
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Nonlinear power spectra, presuming an extreme scale- 
biasing scheme. The density field has been set to 'on' 
;ontaining mass exceeding the largest particle mass in 
Mpc suite, 2.1 X 10^ h~^ A^0, and 'off' everywhere 
: nass cut is most likely lower than anything that could 
) the CfA2 or APM catalogs, except if one assumes an 
small mass-to-light ratio. The result of such a bizarre 
iderjtification scheme is a bias on large scales, due to clear- 
void regions, and an antibias on small scales, due to 
le high peaks in density. We do comparisons with the 
ion suite because the low resolution suite particle mass 



mers, Davis, & Evrard 1995), meaning the power perhaps 
shouldn't be supressed quite as much, and galaxy formation 
will further raise the power. Figure H shows the models' red- 
shift space power spectra, compared to the combined CfA2 



and SSRS2 redshift space power spectrum (ia Costa et al 



1994). Given our choices of model normalization and cosmo- 



logical parameters, the TACDM matter power spectrum is 
nicely consistent with the observed galaxy power spectrum, 
but that leaves no room for galaxy formation or velocity 
bias effects. As for the real-space nonlinear power spectrum 
comparison (figure Vm, this requires significant antibiasing 
for TACDM on scales of 0.3-1 h Mpc"\ Note that under- 
sampling the velocity field will miss the large velocities by 
making the halos physically larger, so it does not make sense 
to perform redshift space comparisons on the large volume 
suite. 

The simplest halo-related quantity to investigate is the 
number density of bound objects eis a function of mass. 
Such 'mass functions' and close relatives such as the X- 
ray temperature function (as in Eke, Cole, & Frenk 199(: , 
for example) are often estimated from the Press-Schechter 
approximation instead of from simulations . Though it has 



beeri checked against s c ale- free simulations ( Efstathiou et al 
198^; Bond et al. 1991; Lacey fc Cole 1994) and against spe- 



cific SCDM ACDM and CHDM mode l s dCarlbcrg fc Couch 



man 1989 Jain fc Bertschinger 1994; Klypin et al. 1995b; 
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Nonlinear redshift-space power spectra 
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Figure 9. Redshift space power spectrum, compared to the com- 



bined C fA2 and SSRS2 redshift space power spectrum (da Costa 
et al. l^gi). Notice that 



1994). Notice that, while ACDM is a good match to this 
power spectrum, there is no room for galaxy formation or velocity 
bias. 



Walter fc Klypin 1996| ; [Bond fc Myers 1996| ), previous stud- 
ies have focused only on a narrow range of masses, typically 
at the cluster scale. With our large simulations, we can check 
the approximation over four orders o f magnitude in mass . 
The Press-Schechter formula we use is Klypin et al. (1995b ), 
equatio ns (1-2). evaluated at z = 0: 



N{> M) 
where 



■ 2 4 f' e[r') 



TtQm Jr a^{r') 



Trexp 



2cr2(r') 



^' (14) 



<r) = 7^ /°° fc^P(fc)iy(fcr) '^^P'"^ dfc, (15) 



27t 



d(fcr) 



(y{r) = —^j k^P{k)W^{kr)dk, 



[sin(2;) — xcos(a;)] tophat 



W{X) = \ fj2/ 



Gaussian 



M 



Qmpcf^O 



(16) 



(17) 



(18) 



and Qm is 47t/3 for a tophat window function, and (27t)''" 
for a Gaussian window function. 

Figure [10 shows the cumulative mass functions esti- 
mated from both suites of simulations, and table H shows the 
Press-Schechter parameters used in that figure. Note that 
in the overlapping region, the two sets of simulation mass 
functions are consistent, and that the high resolution results 
are a significant factor of 1.5-2 below the Press-Schechter 
estimates for all models at the intermediate mass of 10^^ 
h~^ Mq and below. 



Table 3. Press-Schechter fits to simulation mass functions 
Model (5c, t Sc If 



CHDM-2!/ 


1.571 


1.273 


OCDM 


1.693 


1.293 


SCDM b = 1.5 


1.672 


1.236 


TCDM 


1.630 


1.252 


TACDM 


1.732 


1.355 



(5c, t S'lid 5c, g have been chosen to get the same number density 
of clusters with M > 5.5 x 10^'* h^^ Mq as the large-volume 
simulation, for each model. 



This result has been verified recently by other groups. 
Bryan & Norman (1998) see a somewhat stronger discrep- 
ancy at 10^* h~^ Mq, using a spherical overdensity method, 
for cosmological parameters very close to our GHDM-2z^, 
OCDM and SCDM choices (though with s ubstantially larger 



grid c eU sizes for OCDM and SCDM). [SomerviUe et al^ 
(1998) also see an equivalent discrepancy in the differential 



multiplicity function n{M) at z = in the rCDM cosmol- 
ogy, for which the power spectrum is very similar to our 
CHDM-2i^. This result depends upon a completely indepen- 
dent simulation (Jenkins et al. 1997) modeled using adap- 
tive P^M, with halos identified using the Friends-of-Friends 
met hod. The halo mass f unction for the SCDM model from 
the Jenkins et al. (1997 ) simulations is virtually identical 
to ours (G. Lemson, private communication). Given these 
confirmations, we do not believe that the medium-mass dis- 
crepancy we see is an artifact of our simulation method or 
halo finding algorithm. 

As figure O shows, the intermediate and low mass dis- 
crepancy cannot be fixed by adjusting the value of 5c, par- 
ticularly at a mass of ~ 5 x 10^^ h~^ Mq, where th e curves 



cross. Our va lues of 5c, t and (5c, g are consistent with Borgani 
et al. (199"7b| ), except we find that the CHDM-2!/ S^'s are 



not significantly different from the other models. The 5c, t 
values we find for the tophat case are consistent with the 
spherical collapse model. 

The simulation mass functions in figure hcl fall below the 
Press-Schechter predictions for most of their range. For ex- 
ample, in the SCDM high-resolution run, only 60 per cent of 
the particles are within halos with M > 3.4 x 10^"^ h~^ Mq 
&t z = Q. The Press-Schechter prediction is only very slightly 
larger, about 62 per cent for 5c, t = 1.672 and 65 per cent for 
the spherical collapse value (5c, t = 1.686. The mass deficit 
due to smaller abundance of low-mass halos halos in the 
simulations is almost completely compensated for by a small 
excess of very large clusters. The Press-Schechter approxi- 
mation assumes that all the mass must be in halos of some 
size, and this analysis indicates that a significant fraction 
of the mass of the universe should be in small halos. The 
Press-Schechter approximation indicates that for SCDM, as 
much as 20 per cent of the mass is in halos as small as 
10^ h~^ Mq, which is not identifiable in any present-day 
cosmological simulation. The simulations show a significant 
amount of matter that is not in collapsed objects. Most of 
the mass lies in filaments connecting the clusters, many of 
which have only a few identified halos on them. It is con- 
ceivable that much of this mass may be unresolved halos, 
since any A'^-body simulation must have a resolution and/or 
timestep limit below which forces are 'soft,' resulting in dis- 
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Cumulative mass functions 




M (ft-l Mo) 



Figure 10. Cumulative halo mass functions, with Press-Schechter fits. In each panel, the relevant mass functions estimated from the 
two simulation suites is shown by the full curves. Small mass cuts have been applied at M = 3.4 X 10^^ h~^ Mq and 2.2 X 10^^ h~^ Mq, 
for the large and small volume simulations, respectively. Each panel also shows Gaussian (dashed curves) and tophat (dotted curves) 
Press-Schechter mass functions, with Scj. and (5c, g adjusted to agree with the large-volume simulations at 5.5 X 10^* h~^ Mq. The values 
of (5c, t and <5c,g used are given in tabic H. The data points correspond to the observations of BGGMM93 and WEF93, as in figure hi. 



ruption of structure smaller than the limit. Such a mecha- 
nism must be present, since filament halos are necessarily 
not very big, but it is not clear how much of the mass that 
can account for. 

Our two low-no models produce fewer clusters in simu- 
lations than the other models (figure lid). If X-ray temper- 
ature cluster masses are correct, this presents no problem 



for those models. However, if the indications of larger clus- 
ter masses from gravitational lensing are correct, the low- 
f7o models require revision by using less tilt (in the case 
of TACDM) or a larger value of Ho- The former would 
help lessen the disagreement with high-multipole cosmic mi- 
crowave background measurements (figure 0), as a weaker 
tilt would raise the first Doppler peak, but will lead to 
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Press-Schechter mass functions for TACDM 
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Figure 11. Press-Schechter mass functions for TACDM. The high and low density TACDM mass functions from simulations are shown 
in the solid curves. Prom top to bottom at M = 10^^ h~^ A^0i the dashed curves show Press-Schechter mass functions with Gaussian 
filters for (5c, g = 1.0, 1.2, and 1.4, and the dotted curves show tophat filters for (5c, t = 1-4, 1.6, and 1.8. Press-Schechter mass functions 
can be made to agree with our simulations for masses above about 5 X 10^^ h^^ Mq, but not for masses smaller than that. 



the need for even stronger anti-bias to reconcile small-scale 
power with the APM observations. If X-ray temperature 
masses are correct, our parameters for TCDM and CHDM- 
21/ produce too many clusters. For CHDM-2i', the normal- 
ization used here was actually about 10 per cent higher than 
the preferred four-year COBE normalization, so reducing 
the normalization by this factor would probably be enough, 
though this will exacerbate early structure formation prob- 
lems. For TCDM, the only options are to either increase the 



tilt, which is highly disfavored by the small-angle cosmic mi- 
crowave background data as already noted, or further reduce 
the Hubble parameter, which is also strongly disfavored by 
observations. 

The statements above all take the COBE normaliza- 
tion as a fixed constraint. Alternatively, we could turn the 
problem around and use the clusters to determine normal- 
izations and tilts, wit h Hq (and Q.o and O a) as a given. This 
is explored further in Gross et al. (1998). 
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One would now like to investigate statisti cs such as cor- 



rela tion functions, void probability function s (Ghigna et al 
1997 |), p hape statistics (Dave et al. 1997), and other so- 
phisticated statistics. However, to compare to observations, 
we need to know how ma ny galaxies form in each halo. 
Pre vious studies (KNP97 ; [Sfolthenius, Klypin, fc Primack 
1997; Clhigna et al. 1997, for example) have used ad hoc 
'breakup' prescriptions to assign galaxies to halos. We in- 
tend to populate our halos wit h galaxies using a more physi 



Halo mass-weighted autocorrelation functions 



cally m otiv ated approach (as in Kauffmann, Nusser, fc Stein- 
metz 1 )97) based on semi-analytic models including sim- 
plified treatments of gas processes, star f ormation, super- 



nova feedback, and galaxy-g alaxy merging ( Somerville 1997 



Somerville fc Primack 199^). As a result, we do not attempt 



to include any complicated galaxy identification algorithms 
here. 

For certain statistics, one can partially compensate for 
the effect of overmerging by mass weighting. This approach 
is less than ideal because it does not restore the small-scale 
spatial information lost in the overmerging process. Mass 
weighting is equivalent to presuming a halo contains a num- 
ber of galaxies proportional to its mass, and putting all the 
galaxies at the center of the halo. In effect, this clears out 
regions of space around the largest halos' centers, equal to 
their radii, and therefore loses information on scales smaller 
than the largest halo radius (typically 2-3 h~^ Mpc). Since 
very massive halos are rare objects for physically interesting 
cosmoln gical modole, all maeB woightod BtaticticB must bo 



unduly linfliienrpd by sma1l-niimbpr statistical nnis 



We calculate the mass-weighted autocorrelation func- 
tion for the high resolution runs, and the results are shown 
in figure [13. In this figure, a halo mass cut of Af = 3 x 10^^ 
h~^ Mq was used, although the mass weighting makes it 
insensitive to the mass cut. The mass weighting creates a 
spread in the correlation values large enough to prevent 
the test from discriminating among models. To within the 
spread visible in figure h2, all models are roughly c onsistent 
with t he S tromlo-APM autocorrelation function ( Loveday 
et al. 1395). However, there are a few trends visible in the 
figure. SCDM and TCDM are systematically lower in am- 
plitude than the other models, but the effect is not very 
significant given the spread. 



5 CONCLUSIONS 

We have run two suites of simulations with 57 million cold 
particles in boxes of 75 and 300 h^^ Mpc, with the goal 
of studying interesting variants of the CDM family of cos- 
mological models. In this paper, we have made preliminary 
comparisons o f the z = simula tion outputs to data for 
all models. In smith et al. (1998), we used the lower reso- 



lution suite, plus some additional s imulations, to generalize 
the Peacock & Dodds (1994, 1996) procedure for recovery 
of the linear power spectrum corresponding to a given cos- 



molog ical model from observational data. In Wechsler et al 



(1998), we showed that the most massive halos at redshifts 
3, or objects that trace their distribution, can account 



for the obs erved clustering of Lyman-break objects ( Steidel 



et al. 1998 ) for all cosmologies except SCDM. More detailed 
comparisons with observations require assumptions about 
galaxy formation and will be treated in subsequent work. 



CHDM-21' - 

OCDM - 

SCDM 6= 1.5 - 

TCDM - 

TACDM - 

Stromlo-APM f 




r (A.-' Mpc) 

Figure 12. Halo mass-weighted correlation function. Only halos 
with mass above 3 X 10^^ h~^ Mq are included, and all pairs are 
weighted by the product of the two masses. Such a weighting is 
a simple countermeasure for the overmerging problem. The er 



ror ha.rs are the Stromlo-APM autocorrelation function ( Loveday 



et al. 1995) 



Subject to the usual caveats about the uncertainty of galaxy 
formation, we reach the following conclusions in the present 
paper: 

(i) Based on the results of KPH96, who found that ACDM 
models with r2o ~ 0.3 would require strong scale dependent 
anti-bias i n order to be consistent w ith the APM power 
spectrum (Baugh & Efstathiou 1994), we investigated a 
variant of the ACDM model with fio ~ 0.4 and a tilt of 
n = 0.9. We find that this model still requires large antib- 
ias of b^ = PAPMgai/J'dm ~ 0.5 at k — 1 h Mpc~^. Even in 
a simple model in which galaxies are extremely anti-biased 
with respect to dark matter halos, the problem persists on 
scales of r ~ 6h~^ Mpc, because this scale is larger than 
the size of individual halos. To get anti-bias on these scales, 
there would have to be many 'barren' halos containing no 
galaxies. OCDM and TCDM are only shghtly better, stiU 
requiring a strong antibias of fe^ ~ 0.6. Other models con- 
sidered require a weaker antibias at that scale. 

(ii) The TACDM dark matter redshift space power spec- 
trum agrees very well with the redshift space galaxy power 
spectrum from CfA2-fSSRS2 ( [ia Costa et al. 1994| )- This 
leaves no room for the 'positive' bias expected in normal 
galaxy formation, or for velocity biases. For comparison, 
OCDM and TCDM each have room for a modest bias of 
b^ ~ 1.2 at fc = 0.5 h Mpc-\ and CHDM-2!/ and SCDM 
each need b^ ~ 1.5. 

(iii) All models considered here are consi stent with the 
Strom lo-APM real-space correlation function ( Loveday et al. 
1995) on scales of 2-20 h~^ Mpc, largely due to a large 



spread in the model estimates of the correlation function 
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because of mass weighting and small-number statistics for 
large mass objects. 

(iv) The Press-Schechter approximation fits the abun- 
dance of cluster-mass halos very well, with top-hat 
(5c,t=l-57-1.73 and Gaussian 5c, g =1.27-1.35. However, it 
overpredicts the number density of galaxy and small group 
mass objects by a factor of ~ 2, only weakly dependent on 
cosmology, and very weakly dependent on Sc . On mass scales 
of ~ 5 X 10^^ h~^ Mq, it is not possible to compensate for 
the discrepancy by adjusting 5c within reasonable bounds. 

In summary, we conclude that none of the models we 
have investigated can be strongly ruled out by the kind of 
analysis performed here. The CHDM-2i^ model gives the 
best overall agreement with the linear and non-linear tests 
we have considered here, assuming that galaxie s are posi- 
tively bias ed with respect to the dark matter. Gawiser & 
Silk (1998) have shown that a similar CHDM model with 
ill/ = 0.2 in Ni, = 1 neutrino species is a much better fit 
to microwave background and galaxy distribution data than 
any other popular cosmological model. Preliminary analy- 
sis based on the dark matter alone has shown that the re- 
lated CHDM-2z/ model considered in this paper is plausibly 
consisten t with high redshift observations of Lyman-break 
galaxi es ([Wechsler et al. 1990 and damped Lyman-a sys- 
tems ( Klypin et al. 1995aD , but it remains to be seen whether 
this model will produce enough early galaxy formation once 
a more realistic treatment of gas processes and star forma- 
tion is included. More detailed modelling of galaxy forma- 
tion will also be necessary to determine whether the small- 
scale clustering properties of the low-ilo models are indeed 
inconsistent with the observations. In any case we conclude 
that models with flo ~ 0.5 are in better overall agreement 
with the observations t han the lower valu es (fio ~ 0.2-0.3) 
usually considered (e.g. Jenkins et al. 1997). A powerful con- 
straint on Slo, the evolution of cluster abundance with red- 



shift, will be considered in a companion paper (Gross et al 

iggi). 
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